Bulk viscosity, particle spectra and flow 
in heavy-ion collisions 



Kevin Dusling and Thomas Schafer 

Department of Physics, North Carohna State University, 

Raleigh, NC 27695. 

January 12, 2012 

Abstract 

We study the effects of bulk viscosity on px spectra and elliptic flow in heavy ion 
collisions. For this purpose we compute the dissipative correction 5f to the single 
particle distribution functions in leading-log QCD, and in several simplified models. 
We consider, in particular, the relaxation time approximation and a kinetic model for 
the hadron resonance gas. We implement these distribution functions in a hydrody- 
namic simulation of Au-\- Au collisions at RHIC. We find significant corrections due to 
bulk viscosity in hadron px spectra and the differential elliptic flow parameter V2{pt)- 
We observe that bulk viscosity scales as the second power of conformality breaking, 
C ~ i]{c^ — 1/3)^, whereas df scales as the first power. Corrections to the spectra are 
therefore dominated by viscous corrections to the distribution function, and reliable 
bounds on the bulk viscosity require accurate calculations of 5f in the hadronic reso- 
nance phase. Based on viscous hydrodynamic simulations and a simple kinetic model 
of the resonance phase which correctly extrapolates to the kinetic description of a di- 
lute pion gas we conclude that it is difficult to describe the V2 spectra at RHIC unless 
C/s < 0.05 near freeze-out. We also find that effects of the bulk viscosity on the pT 
integrated V2 are small. 



1 Introduction 

One of the fascinating discoveries of the Relativistic Heavy Ion Collider (RHIC) program is 
the near ideal nature of the fluid produced in the collision of two heavy nuclei [1-5]. There 
is a general consensus in the community that the ratio of the shear viscosity to the entropy 
density of the system is no more than a few times the bound t]/s > 1/Att conjectured by 
Kovtun, Son, and Starinets [6]. However, it is difficult to determine the level of accuracy 
that can be obtained when extracting the transport properties. To date, the best estimate 
of the shear viscosity comes from a detailed comparison of particle spectra and elliptic flow 
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Figure 1: Sound speed squared as a function of temperature from the parameterization 
of the lattice QCD equation of state given in [8]. See [9] for a discussion of the various 
parameterizations available for the QCD equation of state. 



with viscous hydrodynamic simulations [7]. But within these state of the art calculations 
there are many systematic uncertainties which are not fully under control. Some of these 
include the precise form of the initial condition, the details of the equation of state, the 
handling of the freeze-out dynamics, and the role of bulk viscosity. Irrespective of its role 
in constraining shear viscosity, the bulk viscosity of the matter produced at RHIC and the 
LHC is clearly an interesting quantity in itself. In this work we will study the effects of bulk 
viscosity on the spectra and the elliptic flow parameter. Our goal is to assess the uncertainty 
in the extraction oi rj/s due to the bulk viscosity, and to identify observables that constrain 
the bulk viscosity. 

The earliest viscous hydrodynamic simulations only included corrections due to shear 
viscosity. One could argue that this may be a safe assumption as there are a number of 
physical systems, possibly relevant to heavy-ion collisions, where the bulk viscosity is zero 
or negligible. For example, it is well known that bulk viscosity vanishes in both the non- 
relativistic and ultra-relativistic limits of a gas when the number of particles are conserved 
[10]. In a weakly coupled quark-gluon plasma, it was found that the bulk viscosity is on 
the order of 1000 times smaller than the shear viscosity [11]. Finally, in the simplest kinetic 
model, the relaxation time approximation, one finds that the bulk viscosity goes as the 
square of the deviation from conformality. 



The above relation was first found by Weinberg for a photon gas coupled to matter [12]. 
It also happens to give parametrically correct results for weakly coupled QCD but not for 
a scalar field theory. In the context of AdS/CFT an analogous relationship [13] has been 
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found, 



C>2r^Q-c^)- (2) 

In this case the bulk viscosity is proportional to the first power of conformal breaking. Based 
on these above examples, it is clear that for a system which is nearly conformally invariant 
(such as weakly coupled QCD) the bulk viscosity will be small. However, lattice QCD 
computations [14] have shown that the equation of state differs strongly from the conformal 
limit at temperatures relevant to heavy-ion collisions (see fig. 1). For example, if the speed 
of sound approaches ^ 0.2 near the phase transition we find ( ~ 0.25?7 using either of the 
expressions (1) or (2) given above. Even larger values ( ~ 0.6r] have been obtained in direct 
lattice studies of the bulk viscosity in the regime T = (1.25 — 1.65)rc [15]. It is therefore 
important to study how bulk viscosity modifies hadronic observables, such as px spectra and 
elliptic fiow. Previous studies of this type can be found in [16-23] . 

We begin by reminding the reader how shear viscosity manifests itself in the spectra of 
produced particles. The equation of hydrodynamics express the conservation of the energy 
momentum tensor, 

d^T^'" = , (3) 

which is given as a sum of ideal and dissipative parts, 

rpt^u = (e + p) uV + Vg^'^ + TTf"" + UAf"" . (4) 

In the above expression for the stress-energy tensor we have used the definition of the three- 
frame projector A^" = g^" + u'^u^ . In the first-order (or Navier-Stokes) approximation the 
dissipative parts of the stress-energy tensor can be written in the local rest frame as 

TT^^' = -T] {d'u^ + cPu' - = -r]a'^ = -2r]{d'u^) , (5) 

n = -Cdku', (6) 

where r] {Q is the shear (bulk) viscosity and (■ ■ ■ ) indicates that the bracketed tensor should 
be symmetrized and made traceless. In principle, it would be satisfactory to solve the 
relativistic Navier-Stokes equations in order to compute the first-order viscous correction 
to particle spectra. However, the first order theory is plagued with difficulties such as 
instabilities and violations of causality. In order to circumvent these difficulties it is necessary 
to use a second order theory, like the one proposed by Israel and Stewart [24,25] or Ottinger 
and Grmela [26,27]. The two theories are qualitatively the same in that they both approach 
the first order theory for small relaxation times. In this work we will not be interested in the 
higher-order corrections arising from the second order theory. Instead we use second order 
hydrodynamics as a practical way to obtain the lowest order correction in going from ideal 
to Navier-Stokes hydrodynamics. 
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The solution to the Navier-Stokes equations will lead to viscous corrections to the re- 
sulting temperature and flow profiles. Particle spectra are then computed using the Cooper- 
Frye [28] formula 

= 77^ [ fiEp)p''da, , (7) 



where is the freeze-out hypersurface taken as a surface of constant energy density in this 
work. For a system out of equilibrium f{Ep) is not the equilibrium distribution function but 
also contains viscous corrections 

f{Ep)^fo{Ep) + Sf{Ep), (8) 

where /o is the usual equilibrium Bose/ Fermi distribution function. The only constraint on 
5f is that the stress-energy tensor remains continuous across the freeze-out hypersurface; 

As shown in [29] this constraint still leaves a lot of freedom in the form of 5/ for shear 
viscosity. It was argued that the functional form of 5/ could fall anywhere between a linearly 
increasing function of momentum to a quadratically increasing function of momentum. These 
two forms of the distribution function lead to qualitatively different behavior for V2{pt) as 
demonstrated by the right plot of fig. 2. By definition V2(pt) is given by 

, , fd(f>cos{2(f)){dN + SdN) 
Jd<l>{dN + 5dN) ' 

where dN is short for dN/ [dpT d(p) and 6dN is the first viscous correction to this. If, as a 
pedagogical exercise we neglect the viscous correction to the distribution function all together 

(which violates energy-momentum conservation across the freeze-out surface), V2{pt) would 
follow the curve labeled '/o' as shown in the left plot of fig. 2. Clearly, the form of the viscous 
correction to the distribution function will play an important role in extracting the shear 
viscosity. 

There is an analogous viscous correction to the distribution function coming from bulk 
viscosity as well. The main goal of this work is to characterize the functional form of Sf 
due to bulk viscosity for various theories and models. We will also show how bulk viscous 
corrections exhibit themselves in spectra as well as some phenomenological consequences. 

2 The Boltzmann transport equation 

Let us first start by setting up the notation that will be used throughout this work. The 
equilibrium distribution functions for bosons and fermions are 
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Figure 2: Typical results for V2{pt) from a viscous hydro dynamic model using a quadratic 
deviation from equilibrium; 5/(p) ~ p'^- The left plot also shows the viscous result without 
including the off-equilibrium correction to the distribution function. The right plot compares 
the quadratic ansatz with a linear ansatz; 5/(p) ~ p- Both curves result in the same shear 
viscosity to entropy ratio. 



where the upper (minus) sign is for bosons and the lower (plus) sign is for fermions. We 
will use capital letters P, Q to label 4-vectors and bold-type p, q for their corresponding 
3-vector components having energy Ep, E^. The magnitude of the three-momentum will be 
written as p, q. The sign convention for the metric tensor is [— , +, +, +] and therefore the 
hydrodynamic fluid four-velocity obeys the normalization condition Ufj,u^ = —1. We also 
use the notation Up = P^(/3)u^(t, x) for the quasi-particle's energy in the laboratory frame 
having four momentum = {P^ = Ep, p) in the local rest frame. 

The starting point for our analysis will always be the Boltzmann transport equation 

Vf{t, X, p) = (dt + Vp-d^ + F- dp) fit, X, p) = -C[/, p] , (12) 

where Vp is the particle's velocity and F is the external force on the particle, 

= dpEp , F = ^ = -d^Ep . (13) 

In this work we will consider only small deviations from local thermal equilibrium and 
therefore expand the Boltzmann equation around the local thermal equilibrium solution 

^> P) = g-/3(..xK(.,x) ^ 1 • (14) 

This procedure is known as the Chapman-Enskog expansion. In the Chapman-Enskog pro- 
cedure we expand the left hand side of the Boltzmann equation in gradients of the thermo- 
dynamic variables and linearize the collision operator in 5/ = / — /eq. Using the following 
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relations^ 



^ - n (l±n (15) 
= np{l±np)P, (16) 



dojp 



the left-hand side of the Boltzmann equation can be written as^ 



np(l ± Up) dj3 



{dt + Vp-d^)/3 + /3{dt + Vp-d^ + F-dp)ujp . (17) 



Let us now assume that the quasi-particles in our system have a dispersion relation of the 
form 

Ep^^/mHP{y:,t)) + p^, (18) 

where we have implicitly included a mass that may be a function of temperature. With this 
dispersion relation the following identities hold 

Making use of the above relations the left-hand side of the Boltzmann equation can be 
rewritten as^ 



/3np{l±np) 2"" ' ^ V3 ' d/3 
where we have defined 

a'^ = 2{d'u^) = (^d'u^ + d^u' - . (21) 

In order to match the kinetic description to hydrodynamics we need to define a covariantly 
conserved energy-momentum tensor in the kinetic theory. There is a subtlety that comes 
about due to the space-time dependence of the mass in the dispersion relation. In order 
to see this, let us first start with the canonical form of the stress-energy tensor which is 
typically used in kinetic theory 

""^^ j J0E~/'^'^^^^' ■ ^^^^ 



^Two useful identities are dup/dp = — np(l ± rip) and d'^Up/dp'^ = np(l ± np)(l ± 2np). 
^Even though we are working in the local rest frame, gradients that are acting on the flow velocity are 
still non vanishing. For example, dfji^ ^ but dfji'^ = since v/^jW,^ = — 1. 

^In deriving this expression we have used the two equilibrium identities dtUi = diln(3 and 9t ln/3 = c^c?jW*. 
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For situations where the dispersion relation is independent of the medium this form is sat- 
isfactory as one can show that energy and momentum is covariantly conserved^ 

di.T^'' = . (25) 

In the case where we have a non-trivial dispersion relation the partial integration can not 
pass through the integration measure. Instead we find that 

Q^J.^,u _ gu ^ ^26) 

where 

We would like to modify the stress-energy tensor such that the above source term vanishes. 
This can be achieved by using the definition 

Throughout this work we will always use this modified form of the stress-energy tensor when 
matching from the kinetic theory to the macroscopic hydrodynamic fields. We stress that if 
the quasi-particle's mass is space-time independent the above two definitions of the stress- 
energy tensor coincide. We also note that these observations are not new. The modified form 
of the stress-energy tensor was used in studies of the bulk viscosity of a hadronic gas [30-33] 
and of scalar field theory [34, 35] . 



3 Relaxation time approximation 

In this section we consider the simplest form of the collision kernel, which is known as the 
relaxation time approximation (RTA) or Bhatnagar-Gross-Krook (BGK) approximation. In 
this model the collision term has the simple form 

C[fM = ^-^^- (29) 



^This can be seen by using the definition of the stress-energy tensor given in eq. (22) and differentiating 
both sides. For the specific case where the dispersion relation is independent of space-time we find 

d^Ti^" = / (^^y^P X, p) . (23) 
In this case the Boltzmann equation is p^dij,f{t,x,p) = —EpC[f,p] and we find 

d,T^" = ^j (^i''^C[/,p]. (24) 
The four-momentum is a colUsional invariant and the right-hand side vanishes. 
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If we define the deviation from equilibrium as Sf{t, x, p) = — /(p) and use the hnearized 
form of the streaming operator given in eq. (20) we find that 



ETT^ ^p(l±^p) 



(30) 



We would now like to identify the relaxation time encoded in 5/ with the transport coeffi- 
cients 7] and C- First we start with the shear viscosity. Looking at any of the off-diagonal 
components of the stress-energy tensor given in eqs. (4) and (28) we find in the local rest 
frame 

5r-^ = -2r^{d^u^) = 1 , (31) 

and the shear viscosity can be identified as 

= 3^ / ^^«(^p)^p(l ± ■ (32) 

If we take a relaxation time of the form^ 

rR(£;p)=ro/3(/3^p)'-" , (33) 
we find the following relation between the shear viscosity and relaxation time 

r 

where the dimensionless phase space integral is worked out in appendix B.l. 

We now come to bulk viscosity, which characterizes the deviation of the pressure from its 
equilibrium value as the fluid expands or contracts more quickly than the time it takes the 
pressure to relax back to its equilibrium value. The bulk viscous pressure, 11, is therefore 
related to the extra pressure from the departure from equilibrium bj. However, the departure 
from equilibrium can not only shift the pressure but also the energy density by an amount 
5e. This shift in energy density will also lead to a shift in pressure, which should not be 
included in the bulk viscous pressure. This is because the bulk viscous pressure should 
only include the difference between the actual pressure and the pressure determined by 
thermodynamics [11] which in our case will be P(e-|-5e). This additional pressure shift must 
therefore be subtracted when defining the bulk viscous pressure^, 

n4r..-.,...)^/^(f-eX^^)./. (3.) 



^We follow the notation of [29] whereby taking a = corresponds to the usual quadratic ansatz. In this 
case the relaxation time grows linearly with momentum, tr ~ -Ep, and X ^ p^. The other extreme case 
follows from a = 1 where now the relaxation is independent of momentum, tr ^ Const., and x ^ v- For 
leading order QCD one numerically finds a = 0.62 and x ^ p^'^^- 

^We have used 

V{eo + Se) « V{eo) + c^Se , (35) 
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E^^^^^ =p' + {m' - ^T') ^p' + m' ^ El . (38) 



Making use of the form of the dispersion relation in eq. (18) it will be convenient to define 
the quantities rh and Ep via 

di(3Ep) , f , dm' 

The following relation between the relaxation time and bulk viscosity coefficient ( then holds, 

C = |Jja(/3m,/3m), (39) 

where the dimensionless phase space integral J^a depends on both the thermal mass m and 
the shifted mass rh. This phase space integral is discussed at length in appendix B.l. In the 
high temperature limit, {T ^ m, m), one finds 



^ = T;^r(6-«), C = i?-rr(6-a) --cM , (40) 



30^^^^'-"^' ^ = ^ 

where the function F, defined in appendix B.l, depends on the statistics of the particles. For 
classical statistics F is the usual Gamma function. From the above formulas we can recover 
the well-known relationship [36] between shear and bulk viscosity, 

C = 15r?Q-c^y . (41) 

We note that this relation is independent of the momentum dependence of the relaxation 
time. 



3.1 Landau matching in the relaxation time approximation 

Landau matching is a way to uniquely specify the energy density e and fiuid four velocity 
in terms of four components of T'*'^. If we use the Landau-Lifshitz convention 

e = u^u,T>'\ (42) 
eu^" = -u^T^''' , (43) 

then the other six independent components of T'^'^ are given by a non-equilibrium stress 
tensor tt^'^ satisfying u^-k^^ = 0. In order that the stress-energy tensor remains continu- 
ous across the freeze-out surface the functional form of 5f must be such that the Landau 
matching condition is satisfied; u^ST^'^ — 0. From eq. (28) the matching condition is 



where from eq. (28) we have 



f d'p ( 



El-T'^]Sf. (36) 
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It is sufficient for the above matching condition to be satisfied in the local rest frame. This 
corresponds to the condition that the shift in energy density stemming from Sf vanishes, 

Let us now look at the energy density shift coming from the off-equilibrium distribution 
given in eq. (30) 

= JJ^> I (ft (I) ' "^(^ * ij - ^-^0 C'^'''"" • ^''^ 

The above expression simplifies considerably when there are no mean fields, Ep ^ Ep, 

S^RTA c I -0^,np{l ± Up) (^^ - clEl^ iPEpf- . (47) 

The above integral vanishes only for a = 1, which is the case where the relaxation time 
TR^Ep) is momentum-independent^. Therefore, if one considers a gas of particles where the 
deviation from conformality comes from the bare mass of the particle only (no mean fields) , 
then the relaxation time approximation can be used if and only if the relaxation time is 
independent of momentum. 

In the presence of mean-fields {i.e. the quasi-particle's mass is temperature dependent) 
we can write eq. (46) as 

iEETA (X J (^«p(l ± »p) (j - '^'K) (f^^p)''" 

In this case taking a = 1 makes the first term vanish, but the second term remains finite 
(even though it may be parametrically small since it is proportional to the coupling). It is 
possible, however, to use the relaxation time approximation consistent with Landau matching 
by a fine-tuning of the parameter a. 



4 Scalar field theory 

The case of a weakly coupled scalar field theory was studied by Jcon [34] where the Boltzmann 
equation and collision kernel were derived from first principles. While the full computation 



'^This is easily seen by using the definition of the sound speed, 
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of the transport coefficients are numerically intensive a lot can be said about the form of the 
off-equilibrium distribution function from certain general considerations. As shown in [35] 
one can compute the transport coefficients in g4>^ + X4>* theory at weak coupling by solving 
Boltzmann equation®, 

(dt + Vp-d^ + F. 9p) fit, X, p) = -C2^2[/, P] - C2^4[/, p] , (50) 

where the coUision operator has been split into a term containing 2 -H- 2 processes and a 
second term involving number changing 2 •<-> 4 processes. While the number changing pro- 
cesses are higher order in the coupling constant (A), they are required in order for a system 
undergoing a uniform expansion or contraction to equilibrate. If number changing processes 
were not included the above Boltzmann equation would have no solution. Formally, this is 
due to the presence of a (spurious) zero mode associated with particle number conservation 
in the 2 ■<-> 2 processes. This zero-mode is not orthogonal to the source term and sub- 
sequently renders the linearized Boltzmann equation non-invertible. We should also point 
out that there is a zero mode corresponding to energy conservation. This zero-mode is not 
problematic since it is orthogonal to the source. 

It is precisely the above behavior of a scalar field theory that allows one to obtain the 
approximate form of the off-equilibrium distribution function. In order to see how this works 
out let us start by linearizing the above Boltzmann equation around its equilibrium solution 

^/(P) = -np{l + np)xn{p)p'p^ {diUj) - np(l + np)Xu{p)dku'' . (51) 

This equation for 6f follows from the Chapman-Enskog expansion eq. 17. The equations in 
the shear and bulk channels can be separated. In the spin (bulk) channel we find 

^p{j- ''^^^^^) = - ' ^^^^ 

where we have written C[Sf, p] to make it explicit that the collision term should be linearized 
around the equilibrium solution. The resulting operators (including the final state symmetry 
factors) are 

C2^2[Sf,p] = h I Tpk^p/k' npnk(l -FnpO(l + nk') 

-'k,p',k' 

X [XM + Xn(^) - Xn(p') - Xn(^')] , (53) 

C2^4[5/, P] = ^ / Tp/k^pk'qq' npnk'nqnq' ( 1 + Up/ ) ( 1 + Uk) 

Jk,p',k',q,q' 

X [Xn(p') + Xn(^) - Xuip) - Xuik') - XuiO) ' Xn(?')] 

- JTjT / Tpk-^p'k'qq' np/nk'^qriq/ (l ^p) (l ^k) 

Jk,p',k',q,q' 

X [Xn (P) + Xn {k) - Xn iv') ' Xn (^') " Xn (?) " Xn (54) 



^For our discussion it will be sufficient to look at a pure A(/)^ theory. 
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where we have used the shorthand J^ — J '^s- The transition rates are given as 



|2 



|2 



- (2^,)(2g,)(2TO,)(2^,)(2g,,) P->'^^'^ + ^ - P' - ^' - g - g').(56) 

Formally, we can solve eq. 52 by inverting the collision operator. Lu and Moore observed 
that the largest contribution will come from the near-zero mode [37] which has the form 

xM=Xo-XiEp, (57) 

where Xi are constants to be determined. Substituting the above form of X-aiv) i'^^'^ ^Pi^ 
channel of the linearized Boltzmann equation, eq. (52), and integrating both sides over all 
phase space we obtain 

Xo = , (58) 

^J- inelastic 



where 



rinelastic = ^pk^p'k'qci' np/nk'^gnq/ (1 + np)(l + 71^} , (59) 

4o Jpkp'k'qq' 



'pkp'k'qq' 

and we have defined the function 



(27r)3Ep V 3 ' op 

which characterizes the deviation of the theory from conformality. The total inelastic cross- 
section given in eq. (59) can be computed by doing the phase space integrals numerically. 
From a phenomenological perspective this is not necessary. Instead, the total inelastic cross- 
section can be related to the bulk viscosity coefficient by using eq. (37). This identification 
leads to 

Xo = |: . (61) 

The constant xi is undetermined by the Boltzmann equation. Instead it is constrained by 
requiring that the deviation from equilibrium does not bring about a shift in the energy 
density. 
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We therefore find the following form for the off-equilibrium distribution function 



XM-Ji^-GE^) , (63) 
where J-" has been defined in eq. (60) and 

= (2^) E^ P pv py ^ 

For completeness, it is worth discussing the parametric behavior of the bulk viscosity at high 
temperature. The bulk viscosity coefficient is given by 

inelastic 

In the high temperature limit we can evaluate T semi-analytically (see appendix B.2). In 
this regime we can ignore the bare and thermal mass of the scalar quasi-particles (up to 
logarithms). The deviation from conformality contained in is controlled by the running 
of the coupling. For a scalar field theory we have 

"^thermal = ^ >m ^ ■ (66) 

3A2 



and using ^(A) = we find that 



16n- 

2rri 



AT ln(7A) l,i5C.(3)/.^ . (67) 

3(327r2)2 '96 ^ ' 

Naively the total inelastic rate would go as A^T^. However, there is a soft enhancement 
which leads to F = ^A^T^ [35]. We therefore find that 

_ A3T3ln^(7A) 

#9(327r2)^ ■ ^^^^ 



5 Leading log treatment in QCD 

In this section we will use the Boltzmann equation in the leading log(T/m£)) approximation. 
In this approximation the dynamics can be summarized by a Fokker-Plank equation which 
describes the momentum diffusion of the quasi-particles. The functional form of Xu ^e 
found by solving a simple ordinary differential equation. We start by discussing the pure 
glue theory and then consider a multi-component QGP. 
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5.1 Pure Glue 



In a leading log approximation, log(T/m£)) is considered to be parametrically large. The 
resulting dynamics describes Coulomb scattering with a small momentum transfer of order 
q ~ gT but with a rapid collision rate of ~ g'^T (up to logarithms). At leading log order 
the linearized Boltzmann equation can be recast as a Fokker-Planck equation [38,39]. This 
equation allows us to determine x(p) in a suitable limit (absence of "gain" terms) by solving 
a differential equation rather than an integral equation. The Fokker-Planck equation is 



Wp(l + '^p) dp 
gain terms 



dp^ 



+ 



np(l + np) 

where is the drag coefficient in the leading log approximation 



dp 



= I^aP , I^A = — ^ log 1 I . 



np(l + np) 



(69) 



(70) 



The Debye mass is given by m|) = \{C a + -^)g'^T'^ with Ca = Nc- Eq. (69) without the gain 
terms is a Fokker-Planck equation for a hard particle undergoing drag and diffusion in a 
thermal bath. In order to conserve energy and momentum the gain terms must be included. 
The gain terms can be written as [38] 



gain terms = 



2^3 



]_d_ 
dp 



p^np{l + rip) 



dE 6 



d_ 
dp 



dP 

Itt 



(71) 



where dE/dt and dP/dt are the energy and momentum transfer to the hard particle from 
the thermal bath per unit time; 



d^ 
'dt 



d^p 
(27r) 



dP 

'dt 



d^p . 

(27r)3J^ 



where 



d 

ip^ -TnAnp{l + np) — 



Sf 



np(l + np) 



(72) 



(73) 



We express the off-equilibrium distribution function in terms of Xn and Xn (^l)- 
Substituting this expression into the Fokker-Planck equation we find that the shear and bulk 
contributions decouple. In the shear shear channel the gain terms vanish and we are left 
with the following ordinary differential equation for Xn{p) 



T 



2np 2\ , 6 



p 



p^ 



(74) 
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Figure 3: Non-equilibrium distribution functions x-n- (red curve labeled shear) and Xn (blue 
curves labeled bulk) of gluons in leading log approximation. The functions x-k and Xn 
defined in eq. (51). We have rescaled x-k by one power of the conformal breaking parameter, 
(1/3-cf), in order to check the expected scaling behavior Xu ~ (l/3-Cs)x7r- The dotted line 
shows the bulk viscous correction Xn before it was made orthogonal to the energy density. 
The curves in this plot were obtained for m^/T = 1, corresponding to a very weak coupling 
as = l/(47r). 



At high momentum (1 + 2np) — )■ 1 and we find 

xAp) = T^P' ■ (75) 

The above differential equation can also be solved numerically. For this purpose two bound- 
ary conditions must be specified. The first boundary condition is that Xw{p = 0) = 0, which 
implies that in QCD soft gluons equilibrate rapidly. The second boundary condition follows 
from the structure of the solution at large momentum. In general the differential equation 
has two independent solutions; one being a polynomial in p and the other growing expo- 
nentially in p. We choose the second boundary condition so that the exponentially growing 
solution is suppressed. In practice, this can be done using a shooting method on x'{p = 0) 
such that x"'{p = Pmax) = 0, which removes the exponential solution. The result of this 
procedure is shown in fig. 3. The shear viscosity can be found using the relation 

^ = Y1 J ^^P*^^ ^ rip)xAp) , (76) 

and we find t]/ (^f^T^ln) = 27.1 in agreement with [39]. 

In the case of bulk (/ = 0) channel, while dP/dt is zero the gain term dE/dt is non- 
vanishing. In order to understand the role of this term we first analyze the Fokker-Planck 
without the gain term 
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Figure 4: Deviation of Quarks and Gluons from equilibrium due to a bulk stress in leading 
log approximation. The dashed curves show the results before the solutions were made 
orthogonal to the energy density. 

This differential equation has one exact zero mode, Xn ^ const., related to particle num- 
ber conservation in 2 o 2 scattering. This zero mode is removed if 2 ^ 3 splitting and 
joining processes are included. We can take this into account by imposing the boundary 
condition XniP = 0) = 0. The second boundary condition is chosen in order to suppress the 
exponentially growing solution as discussed in the shear case. 

The solution obtained in this way is not physically acceptable because it does not respect 
energy conservation. The fact that the collision term conserves energy implies that the 
most general solution of the linearized Boltzmann equation must be of the form Xnip) ~ 
Xn (p) + X^Pj where x^ is a constant and we have used the fact that the leading log collision 
integral is computed using Ep ~ p. It is easy to see that this is a property of the Fokker- 
Planck equation in the bulk channel with the gain term included, but not without it. We 
find that restoring the zero mode Xn P is the dominant effect of the gain term, and that 
Xuip) is very well approximated by the solution of the ordinary differential equation (77). 

The freedom in adding the zero mode has no effect on the calculation of the bulk viscous 
pressure via eq. (37), because any shift in the pressure due to a shift in the energy density 
is projected out. However, in this work we are also interested in the correction to the single 
particle spectra, and in that context the linear term in Xn matters. We therefore fix x^ by 
the requirement that 6f does not contribute to the energy density as required by the Landau 
matching conditions 



In this case there is no need to remove the shift in pressure due to the shift in energy density 
when computing the bulk viscosity 






(79) 
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The numerical solution of eq. (77) is shown in fig. (3). We observe that xb changes sign at 

p ~ 4T, and that for large values of the momentum, p > 7T, the non-equilibrium distribu- 
tion function in the bulk channels scales as the distribution function in the shear channel 
multiplied by one power of the conformal symmetry breaking parameter 

Xn ~ - c^) X. ■ (80) 

Integrating the solution gives (/{T^aD In = 0.44, in agreement with the result in [11]. The 
bulk viscosity scales as the second power of the conformal symmetry parameter, 

C-^ 47.9^-0^^ (81) 

This result has the same structure as the relation obtained in the relaxation time approxi- 
mation, eq. (41), but with a larger numerical coefficient. 



5.2 Quark— Gluon Plasma 

The previous analysis can be easily extended to a multi-component system. For a quark- 
gluon plasma the extension of eq. (77) is [38, 39] 



QA 



p cLa 



2qF{p) 



Closs(x^) + Closs(x^) + - (x'' + X^ - 2x') 
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= Closs(x') - Closs(x^) + - (X^ - X^) 

p 



(82) 
(83) 
(84) 



where x^'^ = Xn^(-P) off-equilibrium distribution functions for gluons and quarks, and 

qi=A,F is the corresponding source term {A adjoint gluons, F fundamental quarks). The 
source and loss terms are different in the shear (/ = 2) and bulk (/ = 0) channels. In the 
bulk channel 



2 \ P 2 ~ 2 



Qiip) = 
Closs(x) = f^iT(^-x"+(^ 
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For comparison, we also show the corresponding source and loss term in the shear channel, 
Qiip) ^ f;, (87) 

/ /I -I- 9r7 9\ fi \ 

Closs(x) = l^iT -x' 
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Figure 5: Differential elliptic ffow of Quarks and Gluons. The solid curves labeled 'Quarks' 
and 'Gluons' represent the quark and gluon elliptic flow using the leading log form of the 
shear viscous correction to the distribution function. In both figures the shear viscosity to 
entropy ratio is rj/s = 0.16. The corresponding dashed curves are the results for a viscous 
hydrodynamic evolution having rj/s = 0.16 and C/s = 0.04. The dashed curves in the left 
plot neglect the bulk viscous correction to the distribution function at freeze-out. The left 
plot should be taken as strictly pedagogical since energy-momentum conservation is violated. 
The right plot shows the complete leading log result. Additional details of the hydrodynamic 
parameters can be found in appendix A. 



The coupled second order differential equations for x^''^ can be solved in the same manner 
as the pure glue case. The result is shown in fig. (4). We observe that there are important 
differences between quarks and gluons, and that there is a shift in the gluon distribution 
due to the presence of quarks. Integrating the distribution functions gives a bulk viscosity 
C/iT^al) In = 0.66 for Nj = 3. 

We are now in a position to compute viscous corrections to the elliptic flow of quarks 
and gluons. Our calculations are based on the 2+1 dimensional second order hydrodynamics 
code described in [2]. See appendix A for details of the hydrodynamic model. We choose an 
initial energy density appropriate for Au + Au collisions at 200 AGeV. The results shown in 
flg. 5 correspond to an impact parameter h = 6.8 fm. The differential elliptic flow parameter 
V2{pt) for quarks and gluons is computed using the strategy outlined in the introduction. 
We have used mo = 2.9T which corresponds to = 0.2. For these parameters leading log 
QCD predicts rj/ s = 0.16 and C/s = 0.08. These values of the transport coefficients lead to 
rather large corrections of the spectra. The results show in flg. 5 were obtained for a smaller 
value of the bulk viscosity, (/^ = 0.04. 

In both the left and the right panel of flg. 5 the elliptic flow parameter V2{pt) in ideal 
hydrodynamics is shown as the solid red hne, and the elliptic flow of quarks and gluons in a 
simulation with shear viscosity only is shown as the solid green and blue curves. The dashed 
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curves in the left panel show the result if bulk viscosity is included in the hydrodynamic 
evolution, but not in the distribution functions (shear viscosity is included in 6f). We note 
that this procedure violates energy-momentum conservation across the freeze-out hyper- 
surface, but it gives an indication of the role that bulk viscosity plays in the hydrodynamic 
evolution. The inclusion of bulk viscosity reduces both the radial flow and the momentum 
anisotropy. These two effects lead to a small reduction of V2{pt) for Pt '^'^ GeV. 

The right panel in fig. 5 shows the full result including the effect of bulk viscosity on the 
distribution function. Comparing with the left panel we clearly observe the importance of 
viscous correction to 5f . From eq. (51) and fig. 4 we can see that the shift in the distribution 
functions due to bulk viscosity is positive at small pr- From fig. 4 the sign change in Xu 
occurs around p/T b. At a decouphng temperature of 150 MeV this corresponds to 
Pt %j 750 McV. Taking into account the boost due to radial expansion the critical px is 
further reduced to pr ^400 MeV, which is barely visible on the plot. At higher momentum 
the bulk viscosity tends to soften the pr spectra. As the spectra enter into the denominator 
in eq. (10) this leads to an increase in V2{pt)- 

Overall, the effect of bulk viscosity on V2{pt) in the regime pr GeV is modest, 
considering that ( is only a factor of four smaller than r). This result is consistent with the 
scaling relations (80) and (81). At very weak coupling ( is suppressed by two powers of the 
small parameter (1/3 — c^), whereas 6f is only suppressed by one power. At strong coupling, 
however, the large numerical coefficient in eq. (81) enhances (/rj relative to Xn/Xn- 

5.3 Leading order behavior at large momentum 

In perturbative QCD the leading order result for the bulk viscosity is governed by small 
angle 2 •<-> 2 scattering, and inelastic 2 •<-> 3 processes are suppressed by a logarithm of 
the couphng constant. At large momenta, p > T/\og{l/g), the logarithmic suppression is 

compensated by the growth of the 2 O 3 reaction with energy. In this regime the correction 
to the distribution function is determined by the physics of energy loss. Arnold et al. showed 
that at leading order in the coupling these effects can included in terms of an effective 1 <H- 2 
collision term [40] 




be 





■(i-x)p) [Xp ~ X^p ~ X(i-j;)p] ) (89) 



be 
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where a,b,c — g,q for quarks/gluons and Xp = Xnip)- The sphtting functions are given 

by 

„ ^ asCAdA I ^ ^2 + n - x)^ ^ + + ~ (90) 

(27r)4 '''^ . . V- (a;(l-a;))V2 
"^^-'^ + + 7^#^ (92) 



7, 



(91) 



7^ 

la 



where k = {2Cf — Ca)/Ca and 

is the transverse diffusion constant that controls energy loss in a quark gluon plasma. We 
can study the effect of the f ^ 2 splitting term on the solution of the Boltzmann equation in 
the bulk channel at large pt- We find that the asymptotic form of Xn is suppressed relative 
to the asymptotic solution for Xi^ by the first power of the conformal symmetry breaking 
parameter, 



Xn = - ^i) XM for p » T\n-\l/g) . 
The asymptotic form of the gluon distribution in the shear channel is given by 



(94) 



Xiip) ^ -^P'/' , (95) 

where we have used Nf — 0. The corresponding result for the quark distribution, as well as 
the dependence on the number of flavors, is given in [29]. 



6 Hadronic Gas 

In the previous section we saw that there are significant differences between the viscous 
corrections to the differential elliptic flow of quarks and gluons. Of course, the spectra of 
quarks and gluons are not directly observable. In this section we study the question whether 
similar differences are expected in the spectra and i'2(pr) of different hadronic species. 

6.1 Low temperature pion gas 

The bulk viscosity of a pion gas was studied by a number of authors [37, 41-43] . Lu and 
Moore argued that the system is similar to the scalar fleld theory studied in section 4, and 
that the bulk viscosity is controlled by number changing processes [37]. We will therefore 
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follow the discussion leading up to eq. (57) and assume that the deviation from equilibrium 
is governed by the near zero-mode, 



X{p) = Xo- XiEp . (96) 

The coefficient xi is determined by Landau matching, and the coefficient xo is controlled by 
the inelastic cross-section, 

XO = , (97) 

inelastic 

where T as written in eq. (60) is a measure of the deviation from conformal behavior. In 
the case of a pion gas wc will ignore mean-field effects (rn-Tr = rn^), and take the deviation 
from conformality to be driven by the bare mass of the pion. In this case J- takes the form 



The total inelastic rate is dominated by the lowest order number changing process which is 
kinematically allowed; tttt TTTTTTTr. The inelastic cross-section also controls the chemical 
equilibration rate of pions. The rate at which a pion chemical potential will return to 
equilibrium is given by [44] 

1 _E,(ferfr,_ jggj 



T-chem. n ' 

where the sum is over all reactions which increase the pion number by Sn^. We can therefore 
make the following identification between the bulk viscosity and chemical relaxation time, 

C = — ■ ■ (100) 

If we use classical statistics, which is valid for m^r 3> T, the phase space integrals appearing 
in J-" can be evaluated analytically. Normalizing the bulk viscosity by the entropy density 
we arrive at the following relationship between the bulk viscosity of a low temperature pion 
gas and the chemical equilibration rate. 



r~ , (101) 



where Ki=i^2,3 is the modified Bessel function of order i — 1,2,3 evaluated at {/3mTr). The 
chemical reaction time arising from inelastic pion reactions can be computed in chiral per- 
turbation theory. For example, the work of [45] (see also [46]) found r^'^'^™- = 450 fm/c at 
T = 140 MeV and r^^'^'^- = 120 fm/c at T = 160 MeV for a pion mass = 138 MeV. Based 
on these calculations we find C/s ~ 0.14 at T = 140 MeV and (/s ^ 0.03 at T = 160 MeV. 
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6.2 Hadronic resonance gas 

The estimate of C for a pure pion gas is likely to be relevant only in a relatively small 
temperature regime. In the regime between the freeze-out and the critical temperature many 
resonances are important. We will assume that the bulk viscosity of a hadronic resonance 
gas is also dominated by number changing processes. If this is the case we may approximate 
the deviation from equilibrium due to bulk viscosity for each hadronic species by the near 
zero-mode 

Srip) = -nlil ± nDdku" (xS - Xi^pJ , (102) 

where Ep^ — v^p^~-i-mf . The coefficient xi (which is the same for all species) is determined 
by the Landau matching condition 

(56 = = / (^^p^^r (P) ' (103) 

where a = n, K, . . . is a sum of all hadronic species in a resonance gas having degeneracy 
i/a- Using the generalization of eq. (37) to a system of multiple species we find 

a 

where 

As in the case of a dilute pion gas we neglect mean-field effects and assume that the deviation 
from conformality is related to the bare masses of the resonances. The off-equilibrium 
distribution in a multi-component system is determined by one parameter, Xi, which is 
common to all species, and A^species parameters Xo ^^^^ different for each species. The 
parameter Xi is determined by the Landau matching condition, and one linear combination 
of the Xo can be related to the bulk viscosity. Exphcit information on inelastic hadronic 
cross-sections is needed to determine the remaining (A'^spedcs — 1) coefficients. 

In this work we will not attempt to compute these inelastic rates. Instead, we will rely on 
a model that is motivated by prior calculations of chemical equilibration rates in a hadronic 
resonance gas [44-48] . Using a phenomenological model for the inelastic cross-section Pratt 
and Haghn showed that the chemical equihbration time near thermal freeze-out is 5 — 10 
times larger for kaons than it is for pions [44]. A similar estimate was also obtained in 
a BUU transport model [47]. Wc therefore expect the bulk viscous correction of kaons to 
be that much larger than pions {i.e. Xo /Xo ~ 5 — 10.). A larger set of resonances (but 
excluding strangeness) was studied by Goity [45] . In this paper the deviation from chemical 
equilibrium (at fixed temperature) is parameterized in terms of effective chemical potentials 
for non-conserved charges like the total number of pions, rho mesons, nucleons plus anti- 
nucleons, etc. Goity finds that the largest relaxation time corresponds to a chemical potential 
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Figure 6: Transverse momentum spectra of pions, protons (left panel), as well as kaons and 
lambda baryons (right panel). The solid lines correspond to shear viscosity only, and the 
dashed lines show the result for shear and bulk viscosity with rj/s = 0.16 and C/s = 0.005. 



for meson (baryon) resonances approximately twice (2.5 times) larger than that of pions near 
the transition temperature. 

In the following we will use the ansatz in eq. (102) and choose Xo each meson and 
baryon species to be a constant multiple Cm and C^, of Xqi 

{Xl Pions 
Cra^Xl Mesons . (106) 
Cfe X Xq Baryons 

Due to the strong p ^ 27c reaction rate we expect the p and vr mesons to be in relative 
chemical equilibrium. This suggests that Pp = 2p^ and therefore Cm ~ 2. Additionally, the 
average pion multiplicity in the strong pp — )■ nvr reaction is n ~ 5 [49], so that 2pi\f ^ 5/i,r 
and therefore Cb ~ 2.5. These numbers are in good agreement with results obtained by 
Goity [45]. The remaining coefficient Xo is related to the bulk viscosity via eq. (104) 

{1 Pions 
Cm Mesons . (107) 
Cb Baryons 

We emphasize that in a complete calculation that includes inelastic rates such as NN — )■ 57r 
the value of ( is completely determined by microscopic dynamics. Without microscopic 
information about inelastic rates we can place bounds on Xo from the observed spectra, and 
then extract bounds on ( from eq. (107). 

Details of the hydrodynamic simulation are described in appendix A. We use the same 
initial conditions and impact parameter as in the case of the pure QGP simulation. The 
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Figure 7: The left panel shows the elliptic flow of pions for a bulk viscosity at freeze- 
out of (C/s)frzout ~ 0.005. The dashed curve shows the result using the linear form of 
the viscous correction given in eq. (102), and the solid curve shows the result using the 
resummed form given in eq. (111). The right panel shows the elliptic flow of pions from 
viscous hydrodynamics when both shear and bulk viscosity are included. The two curves 
labeled 'bulk+shear' are labeled as in the left panel: the dashed line is the linear form of the 
distribution function, and the solid line shows the resummed result. 



equation of state is a parameterization of a lattice QCD equation of state [8]. In the kinetic 
model defined in eq. (102) we include meson/baryon resonances up to a mass of 1.6 GeV 
(mesons) and 1.8 GeV (baryons). We have checked that the corresponding equation of 
state matches the lattice equation of state at freeze-out. Our resonance gas model implies 
xl — — 100C/(sT). We have chosen (C/s)frzout = 0.005, which corresponds to Xp — — 0.5/T. 
Using the average expansion rate {dkU^) at freeze-out the value of Xo can be translated 
into an effective pion chemical at freeze-out, see eq. (109) below. We find yU^r ~ 25 MeV. 
This value is roughly consistent with the pion chemical potential /x^ ~ 10 MeV used in the 
thermal fireball model developed by Rapp [50]. 

We note that we use the same speed of sound, and therefore the same deviation from 
conformality, in our calculations in the quark gluon plasma phase and the hadron resonance 
gas. The difference between the values of C/s in the two phases is connected with the different 
relations between Xu ^"^^ C for the two systems. These relations reflect different physical 
mechanisms for producing bulk viscosity. In the the quark gluon plasma bulk viscosity is 
controlled by momentum rearrangement, and shear and bulk viscosity are intimately related, 
see eq. (81). In the hadron resonance gas model bulk viscosity is dominated by particle 
number changing processes, and there is no direct relation between shear and bulk viscosity. 
The fairly small value of C,/s in the hadron resonance gas is further related to cancellations 
between low-mass and high-mass resonances in eq. (105). 

In fig. 6 we show the px spectra of pions, protons, kaons and lambdas. The shear viscosity 
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Figure 8: Differential elliptic ffow V2{pt) for pions and protons (left panel), as well as kaons 
and lambdas (right panel). The curves are labeled as in fig. 6. The solid lines show the 
result for shear viscosity only, and the dashed lines correspond to shear and bulk viscosity 
with r]/s = 0.16 and C/s = 0.005. 



was chosen to be ///s = 0.16 as in fig. 5. Corrections to the hadronic spectra due to the shear 
viscosity were computed as described in [29]. We observe that, as in the case of quarks and 
gluons, bulk viscosity increases the spectra at small pt, and suppresses the spectra at large 
Pt- The high px suppression is more prominent in the case of pions because the spectra are 
determined by the competition between the constant term Xo and the linear term —xiEp^ 
term, where the constant contribution is bigger in the case of baryons, Xo ^ Xo- The effect 
of bulk viscosity on the elliptic flow parameter V2{pt) is shown in fig. 7. For comparison we 
also show the elliptic flow in the case of an the ideal gas, and in the case of shear viscosity 
only {rj/s = 0.16). We find that bulk viscosity tends to increase elliptic flow for Pt ^ 1 GeV. 
The reason is the same as in fig. 5: bulk viscosity suppresses the single particle spectra at 
large pt, and the spectra enter into the denominator of the definition of V2{pt), see eq. (10). 
The effect becomes very large for p^ > 2.5 GeV. A similar behavior was seen in [17]. Clearly, 
the large pt behavior is unphysical and stems from the fact that the particle distribution 
function becomes negative at some pt- In order to circumvent this we can attempt to do 
a resummation of the viscous correction. We can expand /"(p) to first order in ST and 
chemical potential /i, 

Srip) = <(1 ± nl) (^^ + . (108) 

Comparing this with the form of the off-equilibrium distribution given in eq. (102) we make 
the identification 

/i'^ = -idku')Tx^, (109) 
6T = +{^ku'')T\^ (110) 
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Figure 9: Differential elliptic ffow of pions using the linearized expression for 6f (left) and 
the resummed form of 6f (right). The same V2{pt) can be obtained for Pt ^ 2 GeV when 
increasing rj/s hj a factor of 2.5 as long as the bulk viscosity is increased as well. 



The physics behind this is straightforward. As a system undergoes an expansion (in heavy- 
ion collisions the expansion rate is dkU^ ~ ^) the density of the system drops. However, due 
to the inefficiency of number changing processes there is an excess of particles with respect 
to what would be expected given the energy density of the system. This excess of particles 
can be parameterized by a positive shift in the chemical potential. We can resum the viscous 
correction by using the ideal distribution function with a shifted temperature and chemical 
potential^ 

rip) ^ — ■ (111) 

The above non-equilibrium distribution function is manifestly positive definite. The resulting 
f 2 spectrum is shown in the left panel of fig. 7. At low pt the spectrum matches the linearized 
form, but it has the advantage that it is well-behaved at high pt- 

Resumming the effects of bulk viscosity on the spectra is not as important if shear 
viscosity is also included. Shear viscosity tends to harden the pt spectra, and therefore 
prevents the distribution function from becoming negative (provided 7]/s is sufficiently large). 
In the right panel of fig. 7 we show the elliptic fiow of pions when both shear and bulk viscosity 
are taken into account. In this case we see much better agreement between the linear and 
resummed result even at large pt- We observe that the effect of bulk viscosity on the pion 
V2{pt) is comparable to the analogous correction to the quark V2{pt), despite the smaller 
bulk viscosity used in our simulation of the hadronic phase. This is related to the larger 



^In our calculations we have put the factor in the numerator in order to avoid possible problems 

with Bose condensation in certain regions of phase space. 
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Figure 10: Integrated V2 as a function of the number of participants for pions (left panel) 
and protons (right panel). We show the result in ideal hydrodynamics, the case of only shear 
viscosity with rj/s = 0.16, and the case of both shear and bulk viscosity with rj/s = 0.16 and 
(/s = 0.005. 



numerical coefficient that appears in the relation between ( and (| — c1)x{p) in the quark 
gluon plasma compared to the hadron resonance gas. 

In fig. 8 we compare viscous corrections to the differential elliptic fiow parameter V2{pt) 
for different hadronic species. Reference [29] observed that a simple model for elastic meson 
and baryon cross section reproduces the empirically observed quark number scaling of V2{pt)- 
Fig. 8 shows that bulk viscosity leads to significant modifications of the V2{pt) of individual 
species, but the scaling relations between different species are approximately preserved. 

At a fixed deviation from conformality the off-equilibrium correction to the spectrum 
increases linearly with the bulk viscosity coefficient (. This means that the value of ( cannot 
be increased by very much without resulting in spectra and fiow parameters that are in clear 
disagreement with the data. However, because of the partial cancellation between shear and 
bulk corrections, it is possible to increase both rj and ( simultaneously without changing 
V2{pt) very much. This is demonstrated in fig. 9, where we show that V2{pt %j 2 GeV) is 
fairly stable in the range {7]/s,C/s) = (0.16,0.005) to (r//s,C/s) = (0.4,0.012). 

This result does not imply that the data do not constrain 77 and ( separately. In fig. 10 we 
show the pt integrated fiow parameter f 2 for pions and protons as a function of the number 
of participants. The number of participants was determined from the Glauber model used 
in [2]. We observe that px integrated V2 is quite insensitive to the bulk viscosity. There are 
two reasons for this result. First, for values of C/s in the range studied in this work the effect 
of bulk viscosity on the velocity field is small. Larger values of (/s may lead to stronger 
effects on the integrated V2. Second, because of Landau matching, the pr integrated change 
in the distribution function is small. 
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7 Summary and Outlook 



In this work we examined the functional form of the non-equihbrium correction to the 
particle phase-space distribution caused by bulk viscosity, see the summary in fig. 11. In the 
high temperature quark-gluon phase the distribution function can be computed using the 
leading log approximation. In this limit bulk viscosity is controlled by 2 <H- 2 processes that 
rearrange momentum. Particle number changing 2^3 processes only play an indirect role, 
in that they prevent the development of an effective chemical potential for gluon or quark 
number. 

We showed that there is a significant bulk viscous correction to the quark and gluon 
elliptic fiow even for a fairly small bulk viscosity coefficient. In addition there arc non-trivial 
differences in the quark and gluon off-equilibrium distribution function. These differences 
are related to differences in the transport coefficients and effective masses. While the quark 
and gluon distributions are not directly observable, these distributions serve as direct input 
for calculations of photon and dilepton production from a bulk viscous medium. The effect 
of shear viscous corrections to the distribution function on photon and dilepton production 
was studied in [51-54]. It is conceivable that bulk viscosity is responsible for the large 
elliptic flow of photons as compared to hadrons that was recently observed by the PHENIX 
collaboration [55] . This possibihty is related to the fact that the bulk strain is larger at early 
times, when most photons are produced, and to our observation that bulk viscosity enhances 
V2{Pt) at intermediate pt- 

For the hadron resonance stage near Tc the calculation of the distribution functions is 
more difficult, and one has to rely on simplified models. The simplest model is the relaxation 
time approximation. The relaxation time approximation correctly captures the scaling of 
C and X with the deviation from conformal symmetry, but it cannot predict the functional 
form of xip) (it relates the behavior of x{p) to the unknown energy dependence of r), and 
it is in general not consistent with Landau matching. The relaxation time approximation 
also assumes that shear and bulk viscosity are related to the same process, which need not 
be the case. 

A simple model for theories in which bulk viscosity is controlled by chemical non- 
equilibration is scalar (f)^ theory. In this theory the form of the non- equilibrium distribution 
functions is determined by the exact (energy) and approximate (particle-number) zero modes 
of the collision operator, x — Xo ^ Xi^p- The coefficient of Xo is related to the chemical 
equilibration time r'^'^^™', and xi is fixed by Landau matching. For a given expansion rate 
{d^Uk) — l/r we can also relate x° to the effective chemical potential that describes the 
over-population of the single particle distribution function, ~ — -Xo- 

The bulk viscosity and non-equilibrium distribution function in a low-temperature pion 
gas is correctly captured by the physics of scalar (f)'^ theory with the appropriate chemical 
equilibration time. In this work we assume that this is also true for a hadron resonance 
gas. We assume, in particular, that the non-equilibrium distribution function of the hadron 
species a is of the form x"" — Xo~ XiEpa, where xi is again fixed by Landau matching. The 
relative magnitude of the coefficient Xo ^'^^ different species was fixed by a simple model for 
the effective chemical potentials of meson and baryon resonances. 
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Figure 11: In this figure we summarize different functional forms of tlie correction to tfie 
single particle distribution function due to bulk viscosity, Xn(p)- The curves show the hnear 
and quadratic form of the relaxation time approximation, the result in leading log pure gauge 
theory, and the result in a gas of massive pions. 



In an expanding system inefficiencies in particle number changing processes lead to a 
particle excess, and both xoi^'^Uk) and xiEp{d^Uk) are negative. This means that bulk 
viscosity softens the px spectra of the produced particles. The change in the spectra leads 
to an enhancement of V2{pt) at intermediate momenta ~ (1 — 2) GeV. 

This enhancement tends to cancel against the effects of shear viscosity. We showed, 
however, that the shear viscosity can be determined reliably by focusing on the pt integrated 
elliptic flow parameter. We also showed that bulk viscosity tends to preserve the approximate 
"quark number scaling" observed in in the identified particle V2{pt)- Once r] is fixed bulk 
viscosity is strongly constrained by the spectra and V2{pt)- The main difficulty is that in the 
hadron resonance gas the relationship between x{p) ^i-ud ( is very sensitive to the contribution 
from high lying resonances. 

For the results shown In figs. 6-10 we used (C/'S)frzout ^ 0.005, and found modest bulk 
viscous correction to V2{pt)- In order to obtain a rough bound on the maximum value of (/s 
allowed by the data obtained at RHIC we have studied the dependence of our results on (/s. 
Figure 12 shows the V2{pt) for identified Ks mesons. We have chosen Kg mesons because 
the contribution from resonance decays, which were not included in this work, are negligible. 
Our hydrodynamic model was tuned previously to reproduce the measured spectra using 
shear viscosity only. This implies that the inclusion of bulk viscosity will typically worsen 
the agreement with data. For (C/s)frzout = 0.005 discrepancies with the data are not large, 
and the previous level of agreement could presumably be restored by retuning the parameters 
of the hydrodynamic model. For (C/s)frzout ~ 0.015 the discrepancy with data in the range 
1 ^ Pt ^ 2 GeV is significant, and it is unlikely that agreement with the data could be 
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Pt [GeV] 

Figure 12: Elliptic flow of Ks mesons from viscous hydrodynamics. The hydrodynamic 
model was tuned such that the "shear only" result (solid black curve) fits the data points. 
The short-dashed green curve and long-dashed blue curve show results from viscous hydro- 
dynamics having a bulk viscosity to entropy ratio C/s = 0.005 and C,/ s = 0.015, respectively. 
The data were obtained by the STAR collaboration at RHIC [56] . 



achieved without affecting other observables, like the pt integrated V2. We therefore feel 
that it is safe to claim that the resonance gas model implies (C/s)frzout < 0.015. We plan to 
perform more detailed fits in the future. 

The most important uncertainty in this bound is related to model dependence in the 
relation between x{p) ^"^^ C- the hadron resonance gas this relation depends on the 
inelastic cross-sections of high lying resonances. We can estimate the uncertainty of our 
results by reducing the number of resonances included in the model. For example, if we only 
keep mesons (baryons) with masses below 0.8 (1.0) GeV we find x° ~ — 30C/(sT). This 
relation allows for roughly identical fits to the spectra with a. C,/s larger by about a factor 
of three. We conclude that a more conservative bound is given by (C/s)frzout ^ 0.05. We 
emphasize that the data support a non-vanishing bulk viscosity. Statistical fits to hadronic 
yields [57] show the need to increase the abundance of baryons {i.e. protons + anti-protons) 
through a chemical-abundance factor^° 7^ 1.6 at RHIC energies. This result can be 
naturally accounted for in terms of a non-vanishing bulk viscosity. 

There are a number of issues that we have not addressed in this work. Clearly, more 
work is needed to constrain the bulk viscosity of a hadron resonance gas. We have also not 
taken into account a possible increase in the bulk viscosity near Tc due to critical fluctuations 

^°The abundance factor 7 has to be distinguished from the fugacity A = e''^^ which enhances the abun- 
dance of particles while suppressing that of anti-particles. 
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Figure 13: Bulk viscous pressure —H/C (solid curves) versus proper time along the freeze- 
out hypersurface shown against the Navier-Stokes value dku'' (dashed curve) for (C/'S)frzout ^ 
0.005 (left) and (C/s)frzout ~ 0.015 (right). 



[58,59]. If there is a rapid increase in the bulk viscosity near Tc one also expects a rapid rise 
in the bulk relaxation time. Onuki [60] showed that the bulk relaxation time diverges near 
more rapidly than the bulk viscosity. This implies that the system may free-stream through 
the transition region without significant effects on single particle observables. Clearly, further 
study in this direction is necessary. 
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A Details of the hydrodynamic evolution 

In this appendix we summarize some details of the hydrodynamic calculations that were 
used to compute the velocity and temperature profiles that determine the spectra of pro- 
duced particles. We assume longitudinal boost invariance with initial conditions in the 
transverse plane taken from a Glauber Model (see appendix A in [29] for more details). For 
all non-central collisions we have used an impact parameter of 6 = 6.8 fm, and a decoupling 
temperature Tfeout = 150 MeV. 

We solve second order hydrodynamic equations using a second order fluid model devel- 
oped by Grmela and Ottinger [26,27]. This model is quite similar to the theory of Israel and 
Stewart [24,25]. Grmela and Ottinger introduce a new dynamical tensor variable c^j,. We 
will see below that this variable is closely related to the velocity gradient tensor tt^,^. In the 
local rest frame the stress energy tensor takes the form 

Tl^nF = PiS'' - c^c^') . (112) 
where a is a small parameter, which will be shown to be related to the relaxation time. The 
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Figure 14: Freeze-out hypersurface (Tfeout = 150 MeV) for a central (6 = 0) collision 
with ctq = 0.01 ((C/s)frzout ~ 0.005) shown as the solid black curve and for ctq = 0.03 
((C/'5)frzout ~ 0.015) shown as the dashed blue curve. 



tensor variable c^,^ is conveniently defined to have the property 

c^uu'' = u^. (113) 
We decompose c^y in terms of isotropic and traceless components c and c, 

U^Ui/ -\- -\- C^u , (114) 
Cpiu = ^ (^A - 1) i.9tiu + U^Uy) . (115) 

The equations of motion are dictated by conservation of energy and momentum dpT^^ = 
along with an evolution equation for the tensor variable c^y, 

1 1 

{d\Cpv - dpCxu - dyCpx) = c^i/ , (116) 

In the limit that the relaxation times {tq,T2) are very small the evolution equation yields 

= r2 {d,u^ + d,u' - + '^-t,5'' duu'' . (117) 

Substituting the above equation into T^^pi^p and comparing the result to the Navier-Stokes 
equation the bulk and shear viscosities can be identified as 

7] = T2pa , 

C = lropa. (118) 
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Quarks 



Gluons 



Po 

Co 



2.51 
9.56 X 10-2 
5.25 X 10-1 
8.64 X 10-^ 

1.66 



4.32 
6.28 X 10-2 
9.56 X 10-1 
3.43 X 10-6 

3.48 



Table 1: Parameterization of the leading log QCD off-equilibrium distribution function. We 

use the functional form Xnip) ~ {(^oP^° + Cip^^) In (p/po), for tjid — 3.9 and Nf — 2. The 
above parameterization yields C,/T^ ^ 3.07. 

In our work wc have taken the parameter a = 0.7. These relaxation times arc small enough 
so that the Navier-Stokes limit is approximately maintained near freeze-out. This is demon- 
strated in fig. 13 where the bulk viscous stress 11 is plotted versus the Navier-Stokes expecta- 
tion for a central {b = 0) collision. For reference we also show the corresponding frcczc-out 
hypersurface in fig. 14. The dynamical variable c^i, was initialized to the Navier-Stokes 
value. 

In fig. 5 we show the elliptic flow of quark and gluons obtained in a simulation with a 
pure QGP equation of state. In order to allow for a speed of sound that is different from the 
conformal value = 1/3 we use a polytropic equation of state 



The adiabatic index 7 is chosen in order to fix a constant sound speed = 0.2 compati- 
ble with lattice parametcrizations near T^. The viscous correction to the distribution was 
computed with a Debye mass rriD = 3.9T so that the QGP sound speed is = 0.2, consis- 
tent with the speed of sound used in the hydrodynamic evolution. We employed a simple 
parametrization of the solution of the Fokker-Planck equation for the off-equilibrium distri- 
bution functions. The parametrization is given in table 1. 

All final state hadron spectra shown in this work were calculated using a realistic equation 
of state which is a parameterization of the lattice QCD equation of state from [8]. This 
equation of state matches on to our hadron resonance gas equation of state below T ~ 160 
MeV. The bulk viscosity during the hydrodynamic evolution was assumed to scale with the 
second power of conformality breaking, 



where (Tq is a free parameter chosen to set the desired magnitude of the bulk viscosity 
coefficient near freeze-out. At our freeze-out temperature of 150 MeV the lattice equation 
of state used in this work yields fa 0.15. In section 6.2 we examine a hadronic resonance 
gas with (C/s)frzout ~ 0.005, corresponding to (Tq = 0.01. 



(119) 




(120) 
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B Phase Space Integrals 

B.l Relaxation time approximation 

In the relaxation time approximation wc found the relationship between the shear viscosity 
and energy-dependent relaxation time tr^E^) in eq. (32) which we rewrite here 

/3 



j |^rR(Ep)np(l±np)cip. (121) 



307r2 

If we take a relaxation time of the form 

TR{E^)^T^nPE^f"' , (122) 

the relationship becomes 

V - / (;g^-p(l ± -P) dp . (123) 

Making the change of variables x = /3Ep we find eq. (34) 

V-ro^IaiM (124) 
where the remaining phase space integral is 



oo 



„2 /«_n2n5/2 



X«(/3m) = / ^ i^pJ—n^l ± n^)dx . (125) 

Jl3m ^ 

Even though we have arrived at the above phase space integral by studying the relaxation 
time approximation, it will turn out we will need the same phase space integrals in other 
contexts as well. It is therefore worthwhile to study some limits where analytic results can 
be obtained. For or a classical gas we can replace nx(l ± n^) nx and the phase space 
integral can be computed analytically when a = 

1a=o = lh{l3mfK^{(5m) (126) 

Another case where an analytic expression can be found is in the high temperature limit 
(/3m 0). For a < 4 we find 

X«(/5m = 0) = r(6 - a) (127) 

where for convenience we have defined^ ^ 

r(a;) Maxwell 
T{x) = { r(,T)C+(.T- 1) Bose , (129) 
T[x)C,-{x — 1) Fermi 



^^We have used the relation 

-dx = C±{n)T{n) (128) 



10 T 1 

which can be derived by expanding the numerator in terms of its geometric series and then performing the 
integral of each term in the series individually. The remaining summation will then be of the form 130. 
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for X > 2 and where 



oo /_i_\fc-l 

C±(-)-E4^- (130) 

k=l 

C+(s) is the usual Riemann-Zeta function and C-{^) = (1 ~ 2^~*)C+(s). We will also need 
the a = 4 behavior of of the above phase space integral. For classical and Fermi statistics 
the above results hold as long as we note that linis^i C-(a^) = In 2. We therefore have that 
Xq,=4 = r(2) for classical statistics and Ia=4 = r(2)ln(2) for Fermi statistics. The above 
integral is logarithmically divergent for bosons when a — 4. The divergence is regulated by 
the mass (or thermal mass) of the relevant quasi-particles. We define the following values 
for T{x = 2) 

r(2) Maxwell 
T{x = 2) ^ <( In (f ) - A Bose . (131) 
r(2) ln(2) Fermi 

The relevant phase space integral for bulk viscosity can be found by using the change of 
variable x = /3Ep in eqs. (35) and (36), 

M,m, . £ '^^ - 'fr"'V (l ± nj [1 - c; (l + ^ . (132) 

As in the shear case analytic expressions are available. For a — and classical statistics we 
find 

Ja=o = IsQ-c^^ (/3m)3ir3(/3m) - 6(/3mc,)2 Q - {PmfK^iM 

+ {^mCsY (l3m)Ki(l3m) . (133) 
If both (3m and /3m are taken to zero the integral is 

Ja = Q-^') r(6-«) (134) 

Another limit of interest is when {(3m) — >■ but rh remains finite. This is physically relevant 
since m quantifies the deviations from conformality, which is crucial to keep when studying 
bulk viscosity, while the bare or thermal mass only effects the kinematics in the phase space 
integrals. The only subtlety is if the phase space integral is logarithmically divergent in 
which case the mass serves as a cutoff for the integral. The resulting expression in this limit 
is 

= Q - r(6 -a)- 2{^mcsf Q - ci\ r(4 - a) + {^fhcsf V{2 - a) . (135) 
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B.2 Scalar field theory 



In this section we evaluate the necessary phase space integrals for a scalar field theory. Let 
us first start with the integral labeled T in eq. 60, 



P 



a/3 



n 



p(l + rip) . 



(136) 



In the high temperature limit we can take /3m — )■ while keeping m finite. Using the phase 
space integrals defined in appendix B.l we find^^ 



rp4 

2^ 



(138) 



The function J-" characterizes the deviation from conformality. The relationship between the 
shifted mass m and the sound speed can be found by using the fact that the source term for 
bulk viscosity is orthogonal to the energy-changing zero mode, 







(27r)3 



This leads to 



1 



.2 ~ (fhP)'^^ 



Using this relation we find for bosons 

(mT)2 [l5C+(3) 



67r2 



27r2 



3r(5) 



(139) 
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